package Simplex;
/**
 * ***********************************************************************
 * Compilation: javac Simplex.java Execution: java Simplex
 *
 * Given an M-by-N matrix A, an M-length vector b, and an N-length vector c,
 * solve the LP { max cx : Ax <= b, x >= 0 }. Assumes that b >= 0 so that x = 0
 * is a basic feasible solution.
 *
 * Creates an (M+1)-by-(N+M+1) simplex tableaux with the RHS in column M+N, the
 * objective function in row M, and slack variables in columns M through M+N-1.
 *
 ************************************************************************
 */
public class Simplex {

    private static final double EPSILON = 1.0E-10;
    private double[][] a;   // tableaux
    private int M;          // number of constraints
    private int N;          // number of original variables
    private int[] basis;    // basis[i] = basic variable corresponding to row i
    // only needed to print out solution, not book

    // sets up the simplex tableaux
    public Simplex(double[][] A, double[] b, double[] c) {
        M = b.length;
        N = c.length;
        a = new double[M + 1][N + M + 1];
        for (int i = 0; i < M; i++) {
            for (int j = 0; j < N; j++) {
                a[i][j] = A[i][j];
            }
        }
        for (int i = 0; i < M; i++) {
            a[i][N + i] = 1.0;
        }
        for (int j = 0; j < N; j++) {
            a[M][j] = c[j];
        }
        for (int i = 0; i < M; i++) {
            a[i][M + N] = b[i];
        }

        basis = new int[M];
        for (int i = 0; i < M; i++) {
            basis[i] = N + i;
        }

        solve();

        // check optimality conditions
        assert check(A, b, c);
    }

    // run simplex algorithm starting from initial BFS
    private void solve() {
        while (true) {

            // find entering column q
            int q = bland();
            if (q == -1) {
                break;  // optimal
            }
            // find leaving row p
            int p = minRatioRule(q);
            if (p == -1) {
                throw new RuntimeException("Linear program is unbounded");
            }

            // pivot
            pivot(p, q);

            // update basis
            basis[p] = q;
        }
    }

    // lowest index of a non-basic column with a positive cost
    private int bland() {
        for (int j = 0; j < M + N; j++) {
            if (a[M][j] > 0) {
                return j;
            }
        }
        return -1;  // optimal
    }

    // index of a non-basic column with most positive cost
    private int dantzig() {
        int q = 0;
        for (int j = 1; j < M + N; j++) {
            if (a[M][j] > a[M][q]) {
                q = j;
            }
        }

        if (a[M][q] <= 0) {
            return -1;  // optimal
        } else {
            return q;
        }
    }

    // find row p using min ratio rule (-1 if no such row)
    private int minRatioRule(int q) {
        int p = -1;
        for (int i = 0; i < M; i++) {
            if (a[i][q] <= 0) {
                continue;
            } else if (p == -1) {
                p = i;
            } else if ((a[i][M + N] / a[i][q]) < (a[p][M + N] / a[p][q])) {
                p = i;
            }
        }
        return p;
    }

    // pivot on entry (p, q) using Gauss-Jordan elimination
    private void pivot(int p, int q) {

        // everything but row p and column q
        for (int i = 0; i <= M; i++) {
            for (int j = 0; j <= M + N; j++) {
                if (i != p && j != q) {
                    a[i][j] -= a[p][j] * a[i][q] / a[p][q];
                }
            }
        }

        // zero out column q
        for (int i = 0; i <= M; i++) {
            if (i != p) {
                a[i][q] = 0.0;
            }
        }

        // scale row p
        for (int j = 0; j <= M + N; j++) {
            if (j != q) {
                a[p][j] /= a[p][q];
            }
        }
        a[p][q] = 1.0;
    }

    // return optimal objective value
    public double value() {
        return -a[M][M + N];
    }

    // return primal solution vector
    public double[] primal() {
        double[] x = new double[N];
        for (int i = 0; i < M; i++) {
            if (basis[i] < N) {
                x[basis[i]] = a[i][M + N];
            }
        }
        return x;
    }

    // return dual solution vector
    public double[] dual() {
        double[] y = new double[M];
        for (int i = 0; i < M; i++) {
            y[i] = -a[M][N + i];
        }
        return y;
    }

    // is the solution primal feasible?
    private boolean isPrimalFeasible(double[][] A, double[] b) {
        double[] x = primal();

        // check that x >= 0
        for (int j = 0; j < x.length; j++) {
            if (x[j] < 0.0) {
                System.out.print("x[" + j + "] = " + x[j] + " is negative");
                return false;
            }
        }

        // check that Ax <= b
        for (int i = 0; i < M; i++) {
            double sum = 0.0;
            for (int j = 0; j < N; j++) {
                sum += A[i][j] * x[j];
            }
            if (sum > b[i] + EPSILON) {
                System.out.print("not primal feasible");
                System.out.print("b[" + i + "] = " + b[i] + ", sum = " + sum);
                return false;
            }
        }
        return true;
    }

    // is the solution dual feasible?
    private boolean isDualFeasible(double[][] A, double[] c) {
        double[] y = dual();

        // check that y >= 0
        for (int i = 0; i < y.length; i++) {
            if (y[i] < 0.0) {
                System.out.print("y[" + i + "] = " + y[i] + " is negative");
                return false;
            }
        }

        // check that yA >= c
        for (int j = 0; j < N; j++) {
            double sum = 0.0;
            for (int i = 0; i < M; i++) {
                sum += A[i][j] * y[i];
            }
            if (sum < c[j] - EPSILON) {
                System.out.print("not dual feasible");
                System.out.print("c[" + j + "] = " + c[j] + ", sum = " + sum);
                return false;
            }
        }
        return true;
    }

    // check that optimal value = cx = yb
    private boolean isOptimal(double[] b, double[] c) {
        double[] x = primal();
        double[] y = dual();
        double value = value();

        // check that value = cx = yb
        double value1 = 0.0;
        for (int j = 0; j < x.length; j++) {
            value1 += c[j] * x[j];
        }
        double value2 = 0.0;
        for (int i = 0; i < y.length; i++) {
            value2 += y[i] * b[i];
        }
        if (Math.abs(value - value1) > EPSILON || Math.abs(value - value2) > EPSILON) {
            System.out.print("value = " + value + ", cx = " + value1 + ", yb = " + value2);
            return false;
        }

        return true;
    }

    private boolean check(double[][] A, double[] b, double[] c) {
        return isPrimalFeasible(A, b) && isDualFeasible(A, c) && isOptimal(b, c);
    }

    // print tableaux
    public void show() {
        System.out.print("M = " + M);
        System.out.print("N = " + N);
        for (int i = 0; i <= M; i++) {
            for (int j = 0; j <= M + N; j++) {
                System.out.printf("%7.2f ", a[i][j]);
            }
            System.out.print("\n");
        }
        System.out.print("\n value = " + value());
        for (int i = 0; i < M; i++) {
            if (basis[i] < N) {
                System.out.print("\nx_" + basis[i] + " = " + a[i][M + N]);
            }
        }
        System.out.print("\n");
    }

    public static void test(double[][] A, double[] b, double[] c) {
        Simplex lp = new Simplex(A, b, c);
        System.out.print("value = " + lp.value());
        lp.show();
        double[] x = lp.primal();

        for (int i = 0; i < x.length; i++) {
            System.out.print("x[" + i + "] = " + x[i] + "\n");
        }
        double[] y = lp.dual();
        for (int j = 0; j < y.length; j++) {
            System.out.print("y[" + j + "] = " + y[j] + "\n");
        }
    }

    public static void test1() {
        double[][] A = {            
            {1, 1},
            {2.5, 11},};
        double[] c = {0.8, 1.3};
        double[] b = {130,250};
        test(A, b, c);
    }
}
